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HIGH-RESOLUTION  ARRAY  WITH  PRONY,  MUSIC,  AND 

ESPRIT  ALGORITHMS 


1.  INTRODUCTION 

This  report  investigates  the  array  high-resolution  properties  of  three  algorithms:  Prony, 
MUSIC  (Multiple  Signal  Classification),  and  ESPRIT  (Estimation  of  Signal  Parameters 
via  Rotational  Invariance  Techniques). 

The  Prony  algorithm  [1]  for  estimating  the  parameters  of  sinusoids  in  noise  is  very 
old.*  (Its  discoverer,  the  Baron  de  Prony,  lived  around  the  time  of  the  French  revolution.) 
Unlike  MUSIC  and  most  other  eigenspace-based  methods,  it  can  work  with  correlated 
signals,  but  because  of  the  assumed  signal  model  (sinusoids  in  noise),  it  is  geometry- 
dependent.  Specifically,  it  must  assume  that  the  array  is  linear,  and,  therefore,  the  range  of 
applicability  is  greatly  reduced.  Nevertheless,  the  importance  of  resolving  correlated 
signals  makes  it  an  interesting  high-resolution  algorithm. 

The  MUSIC  [2]  algorithm  has  received  a  lot  of  attention  due  to  its  excellent  behavior 
when  the  signals  have  even  a  small  degree  of  decorrelation.  However,  it  breaks  down 
when  the  signals  are  perfectly  or  nearly  perfectly  correlated  (in  the  statistical  sense 
explained  below).  Unfortunately,  this  can  occur  in  many  cases  of  interest,  especially  with 
multipath,  and  it  is  then  difficult  to  distinguish  multipath  from  true  targets. 

ESPRIT  [3]  is  a  recently  developed  eigenspace-based  technique  that  has  the  same 
excellent  resolution  properties  as  MUSIC,  but  is  much  more  computationally  efficient.  It 
can  also  be  adapted  to  the  correlated  case. 

This  report  concentrates  primarily  on  an  investigation  of  the  Prony  algorithm  in  the  case 
of  perfect  correlation.  We  also  present  many  results  from  MUSIC  in  the  uncorrelated  and 
nearly  correlated  cases  for  comparison.  In  Appendix  A,  we  take  a  brief  look  at  ESPRIT  for 
the  perfectly  uncorrelated  case.  In  Appendix  B,  we  present  some  results  from  a  variation 
of  MUSIC  called  "spatially  averaged"  MUSIC  that  can  distinguish  between  correlated 
signals.  Indeed,  the  results  are  comparable  to  Prony  in  the  case  of  perfect  correlation. 
Finally,  in  Appendix  C,  we  glance  at  some  of  the  recent  literature  on  MUSIC. 

2.  THEORY 

We  will  now  describe  the  theory  behind  the  first  two  algorithms  examined  in  this 
report:  the  Prony  algorithm  [1]  and  the  MUSIC  algorithm  [1,  2,  4j.  In  all  cases  we 
assume  that  the  receiving  antenna  is  in  the  far  field  of  a  limited  number  of  discrete 
narrowband  sources,  so  that  at  any  given  instant,  the  incoming  signal  is  closely 
approximated  by  a  set  of  discrete  plane  waves  coming  from  different  directions  of  arrival 
(DOA).  (The  definition  of  narrowband  in  this  case  is  that  the  time  of  travel  of  a  wavefront 
across  the  array  is  much  less  than  the  duration  of  the  pulse.)  We  consider  the  one¬ 
dimensional  case,  i.e.,  the  array  is  one-dimensional  and  the  angle  of  arrival  varies  in  the 
horizontal  plane  only,  measured  with  respect  to  boresight.  For  simplicity,  we  ignore 
coupling  among  the  elements,  and  the  elements  are  assumed  to  all  be  the  same. 
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Prony  s  original  algorithm  was  an  exact  fit  to  the  noiseless  case.  With  overdetermined  equations,  we  can  make  i 
least-squares  fit  to  the  noisy  case  This  is  equivalent  to  the  popular  Modified  Covariance  algorithm  of  spectral 
estimation  [  1 ,  p.  314], 
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At  any  given  instant,  there  is  a  set  of  complex  baseband  outputs  across  the  elements, 
which  will  be  referred  to  as  a  snapshot.  A  set  of  discrete  snapshots  is  recorded  over  some 
relatively  short  interval  of  time  during  which  the  number  of  signals  and  their  angles  of 
arrival  are  assumed  to  remain  essentially  constant.  Also,  the  statistics  of  the  received 
signals  must  remain  stationary  during  this  interval.  This  includes  the  correlations  between 
the  (complex)  amplitudes  of  the  different  plane  waves  (the  different  signals  may  have 
randomly  varying  amplitudes  and/or  initial  phases  from  snapshot  to  snapshot),  the 
correlations  between  the  receiver  noise  at  the  elements  (each  element  has  its  own  receiver 
and  the  antenna  is  assumed  to  be  internally  noise  limited),  and  the  correlations  between  the 
(complex)  amplitudes  and  the  noise  at  the  elements.  (See  Ref.  5  for  further  explanation.) 

2.1  Prony  Algorithm 

The  Prony  algorithm  models  the  signal  across  the  elements  as  a  sum  of  complex 
exponentials.  This  model  is  correct  only  when  the  array  is  linear,  and  therefore  the  Prony 
algorithm  is  geometry-dependent  (unlike  MUSIC  and  other  eigenspace-based  algorithms). 
In  this  respect,  it  is  inferior  to  MUSIC;  however,  it  has  the  advantage  over  MUSIC  of 
being  able  to  resolve  perfectly  correlated  signals  (which  may  occur,  for  example,  in  the 
case  of  multipath).  Also,  the  Prony  algorithm  degenerates  when  noise  is  added,  as  this 
again  implies  a  violation  of  the  assumed  signal  model.  This  signal,  for  a  given  snapshot,  is 
represented  as 


sn(tj)=£  cm(tj)ein(kDsin9m)  n  =  0,1,2  •••  (1) 

m  =  l 

where  tj  is  the  time  of  the  j1*1  snapshot,  sn  is  the  value  at  the  n1*  element,  K  is  the  number  of 
incoming  signals  (narrowband  plane  waves),  cm  is  the  instantaneous  amplitude  of  the  m1*1 
signal  at  the  0th  element  (the  element  pattern  is  lumped  into  this  constant),  n  is  the  element 
number,  k  is  the  wavenumber,  D  is  the  spacing  between  elements  (A./2  in  all  cases 
considered  in  this  report),  and  0m  is  the  angle  of  arrival  of  the  mlh  signal. 


K 


With  the  following  compact  notation 

vm  =  kD  sin  0m 


—  [Ci  C2  •  •  •  Ck]^ 
dj  =  [e'v‘n  eivzn  •  •  •  eiVKn]T 


(2) 


we  can  rewrite  the  first  equation  as 


rv 

Sn  =  CT  dn  =  X  (fj) ein  (kD  sin  0m)  =  X  (tj)  einVn 


m-1 


m=  1 


By  introducing  a  shifting  operator 


e1V1  0 


D  = 


0 


e'V2 


we  can  write 


0  0 


0 

0 

e'Vk 


(3) 


(4) 
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Sn  —  CT  dj] 

Sn-l  =  £T  dn-1  =  CT  D1  dn 
s„.j  =  £TDjdn 


which  can  be  written  more  concisely  as 

Xl  =  [sn  sn-l  sn-2  •  • '  Sn-(L-1)] 


with 


and 


\  - 


Cl 

C2 

CK 

— iv. 

— i  VS 

C, 

e  1 

c2  e  2 

••  cKe  ‘ 

-i2v, 

-iiv. 

— \2 

Cl 

e  1 

c2  e  2 

■■  cKe 

Cl 

e-i(L-l)v, 

c2e-(L-,)^  • 

-  cK  e'i(L_1)V| 

(5) 


(6) 

(7) 


(B) 


(L  is  some  integer  such  that  L-l  <  n  <  N-l  and  K<L.  Its  purpose  will  become  clearer 
below.) 


Now  any  rectangular  matrix  can  be  given  a  singular  value  decomposition  (SVD)  [1], 
By  applying  this  to  the  matrix  we  can  eventually  arrive  at  an  algorithm  to  determine  the 
directions  of  arrival  0m,  as  follows: 


Ct 

£tD  1 


y  v 
zk  v 


H 

K 


CT  D  L+1 


Uk  Uk  =  Ikk 

UkUk*1kk  (9) 

VK  =  Ikk 

where  U  is  an  L  by  K  matrix,  Z  is  a  K  by  K  diagonal  matrix,  and  V  is  a  K  by  K  matrix. 
The  diagonal  elements  of  Z  are  the  singular  values. 

We  can  then  use  linear  algebra  to  determine  dn  in  terms  of  Xn 
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Xn  =  (Uk  2k  Vk)  dn 
(substitution  for  £) 

UgXn  =  (U^Uk)Zk  dn  =  lKK2KV^dn 
(premultiply  by  U{<) 

VkIJc1  U^XJ1  =  (VKV»)dJ1  =  £ in 

(premultiply  by  Vk  2k) 

which  gives 

dn  =  V1Xn 

where  is  the  Moore-Penrose  pseudoinverse  defined  by 


(10) 


(11) 


£  1  =  Vk  2k  Uk  •  (12) 

Using  this  expression,  we  can  then  determine  a  one-step  predictor,  to  express  Xn+j  in 
terms  of  Xn 


Xn+j  =  ^  D  V1  X„  =  A  X* 

(13) 

(substitution  into  Xn+i  =  ^  dn+i  =  £  Dd„) 

where  A  =  £  D  is  called  the  system  matrix. 

The  system  matrix  satisfies  the  following  eigenvalue  equation 

A£i  =  ^ 

X.j  =  e'vi  i  =  1,2,  ...  K  (14) 

Xi  =  0  i  =  K+l, ...  L 

as  can  be  determined  by  simple  substitution  for  A  and  the  fact  that  the  pseudoinverse  is  a 
left  (though  not  right)  inverse  to  itself: 


A  £  =  £  D  (4'1  5)  =  5D- 

(15) 

(remember  that  D  is  diagonal) 

(Notice  that  the  directions  of  arrival  are  found  from  the  eigenvalues;  they  are  the  arcsines  of 
Vi/kD  for  i  =  1....K.) 
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We  can  represent  the  system  matrix  as  follows 

2£n+l  =  A  Xn  (L-l  <  n  <  N-l) 


Sn+1 

‘A1 

A2 

a3  . 

•  al 

Sn 

S 

n 

1 

0 

0  . 

.  .  0 

s  , 

n-l 

S  . 
n-l 

0 

1 

0  . 

.  .  0 

S  0 
n-2 

• 

0 

0 

1  . 

.  .  0 

.  ‘  j 

0 

0 

0  . 

1  0 

L  L 

Sn  +  1  —  Sn-(k-l)  Sm+L  =  ^  Sm  +  L-k  HI  =  0, 1, 2, ...  N  -  L  -  1  .  (16) 

k  =  l  lt  =  ] 

( fhe  form  of  the  matrix  is  clear  from  the  following  considerations:  The  first  row  of  the 
matrix  consists  of  the  coefficients  that  give  Sn+1  as  a  linear  combination  of  the  previous  L 
s-values,  and  the  remaining  rows  of  the  matrix  are  simply  identities  among  the  Sn.) 

The  Ak  can  be  solved  in  a  least  squares  sense  by  the  normal  equations  of  linear  algebra. 
Specifically,  the  last  set  of  equations  above  has  the  form 


(with  y  corresponding  to  the  Ak)  v 

and  an  optimal  (in  the  least-squares  sense)  y  can  be  found  by  minimizing  the  functional  x  - 
C  x  w^h  respect  to  y 


min  [x  -  C  yj2  =>  CHx  =  CHCy, 


(18) 


which  can  be  re-expressed  as 

N-L-l  *  L  N-L- 1  * 

X  sm+L-j  —  X  X  (sm+L-i  ^m+L-k)  Ak  •  (19) 

m  =  0  k-1  m=0  J 

By  defining  the  correlation  matrix 


R  =  I  X^Xj 

n=L 

IT  =  [Roo  Rio  •  •  ■  Rn-L-1,  o] 


(20) 


and  the  vector 


A  =  [A]  A2  Aj  ■  •  •  Al] 


(21) 


we  can  rewrite  the  above  equations  concisely  as 

RjO  =  X  Rjk 

k  =  l 


(22) 
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or 


I  =  R  A. 


(23) 


We  can  then  solve  for  the  vector  A 


A  =  R  1 1 

(provided,  of  course,  that  R  is  non-singular) 


(24) 


which  readily  gives  the  system  matrix  A.  By  solving  the  eigenvalue  equation  of  the 
system  matrix,  we  can  finally  determine  the  directions  of  arrival  from  the  eigenvalues,  as 
described  above. 


An  alternative  solution  that  lessens  the  computational  load  by  projecting  the  problem 
onto  a  smaller  subspace  is  as  follows.  The  correlation  matrix  has  an  SVD  decomposition 

R  =  VK  £k  Vg  (25) 

where  the  right-hand  matrix  is  the  Hermitian  of  the  left-hand  matrix  because  R  itself  is 
Hermitian.  We  then  project  the  data  vectors  onto  the  K-dimensional  subspace  spanned  by 
the  singular  vectors  that  comprise  the  columns  of  Vk 


Xn=  VjUn 

which  gives  a  reduced  correlation  matrix 


(26) 


R  =  I  x;xl=  I  v£x;xJvk  =  v»rvk  =  2:k 


N-l 


n=L 


n=L 


(the  last  step  results  from  various  forms  of  the  relation  Vk  =  1  kk) 


(27) 


equal  to  the  matrix  of  singular  values 

R  =  Ik. 

We  can  then  determine  a  reduced  one-step  prediction  (or  system)  matrix: 

Xn+l  =  A  Xn 
(given  previously) 

VkXh^AVkX* 

(from  the  definition  of  the  reduced  data  vector) 

vJVKXn+1=V{  A  VkXh 
(premultiplying  by  VK) 

Xn+!  =  AX* 

a=v£aVk. 

(and  using  v£vr*  =  1KK) 


(28) 

(29) 

(30) 

(31) 

(32) 
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Now  defining  the  matrix  Ri 


N-2  — T  N-2  — *  ~T  — *  ~  — *  _ 

Rl=  I  Xn+l  Xn  =  I  A  XnXrl=A  R  =  A  IK 

n=L- 1  n=L-l 

(every  step  a  substitution  of  a  previous  identity) 


(33) 


we  can  express  the  reduced  system  matrix  as 


A  - R  i  Ik  . 

(taking  the  conjugate  of  both  sides  and  using  the  diagonali^y  of  Ik) 


(34) 


It  can  be  shown  that  the  reduced  system  matrix  above  has  the  same  nonzero  eigenvalues 
as  the  full  system  matrix.  However,  its  dimension  is  only  K  by  K,  so  that  to  find  the 
eigenvalues,  we  need  solve  only  a  K-dimensional  characteristic  polynomial  instead  of  a 
possibly  much  larger  one.  (The  eigenvalues,  we  remember,  give  the  directions  of  arrival.) 


a 

^■i  eivi  (i  =  I,...,  K)  S  =  Vj[£i. 


(35) 


In  the  case  of  noise,  the  Prony  method  degenerates,  as  the  assumed  sinusoidal  model 
only  approximates  the  data.  We  now  need  some  kind  of  averaging,  and  there  are  different 
ways  of  doing  this.  (We  assume  that  the  noise  is  zero-mean  and  stationary.) 

The  Prony  space  method  calculates  the  directions  of  arrival  (from  the  eigenvalues) 
separately  for  each  snapshot  and  then  averages  them  over  all  the  snapshots.  Because  of  the 
assumed  stationarity  and  zero  mean  of  the  noise,  the  averaged  estimates  of  the  directions  of 
arrival  should  asymptotically  approach  the  true  values.  In  averaging  the  DOAs,  we  must 
assume  that  the  rank  K  (the  number  of  singular  values)  remains  constant  from  snapshot  to 
snapshot,  or  else  we  must  make  some  provision  to  account  for  changing  rank.  Also,  the 
variance  of  the  eigenvalues  cannot  be  so  great  that  we  cannot  identify  related  eigenvalues 
from  snapshot  to  snapshot.  However,  if  the  noise  is  this  great,  the  method  probably  would 
not  work  anyway. 


The  Prony  space-time  method,  on  the  other  hand,  first  averages  the  relevant  matrices 
over  all  snapshots  and  then  determines  the  eigenvalues  from  the  averaged  matrices.  There 
are  different  ways  of  doing  this.  At  first,  I  averaged  R  and  Ri,  performed  an  SVD  on  the 
averaged  R,  calculated  the  inverse  of  I  from  these  singular  values,  and  then  used  this 
inverse  along  with  the  average  of  Ri  in  the  equation 


A  =  Ri  Ik  (36) 

to  find  the  reduced  system  matrix.  Then  I  calculated  the  eigenvalues  from  this  reduced 
system  matrix. 

This  did  not  work  very  well.  (It  gave  poorer  results  than  the  Prony  space  algorithm, 
and  averaging  did  not  seem  to  help  much.)  I  then  tried  an  alternative  technique,  which 
seems  to  work  somewhat  better  than  the  Prony  space  algorithm,  and  the  results  of  which 
are  presented  in  this  report.  This  consisted  simply  of  averaging  R]  and  the  inverse  of  I 
calculated  from  each  snapshot,  and  then  using  these  two  averaged  matrices  in  the  equation 
above.  (I  did  not  use  the  average  of  the  correlation  matrix  R.) 
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2.2  MUSIC  Algorithm 

The  MUSIC  algorithm  has  two  advantages  over  the  Prony  algorithm.  It  is  geometry- 
independent  (i.e.,  it  does  not  require  the  array  to  have  a  specific  geometry)  and  it 
incorporates  noise  into  its  mathematical  structure  to  reduce  its  effects.  However,  unlike 
Prony,  it  breaks  down  in  the  case  of  perfectly  or  nearly  perfectly  correlated  signals,  and 
this  case  is  often  of  interest  (for  example,  in  distinguishing  multipath  from  true  targets). 
(The  precise  definition  of  "correlated  signals"  will  become  clear  in  the  discussion.) 

MUSIC  starts  by  explicitly  incorporating  noise  into  its  signal  model 

z„  (tj)  =  £  Wn  (0m)  cm(tj)  +  Vn(tj)  n  =  0,  1 . N- 1  (37) 

m-1 

The  various  quantities  are  as  follows:  zn  is  the  signal  at  element  n  at  time  tj.  cm  is  the 
complex  amplitude  of  plane  wave  m  at  reference  element  0  at  time  tj.  Wn(0m)  is  the 
steering  vector  (at  element  n  in  direction  0m)  that  accounts  for  both  the  delays  of  the 
various  plane  waves  between  the  reference  element  and  element  n  and  also  any  weights  that 
may  be  applied  to  the  elements.  (It  is  the  explicit  incorporation  of  this  steering  vector  that 
allows  independence  from  geometry.)  Finally,  vn  is  the  noise  at  element  n  at  time  tj.  It  is 
assumed  that  the  noise  is  uncorrelated  from  element  to  element  and  is  also  uncorrelated  with 
the  signal. 

The  signal  model  can  be  written  more  concisely  as 

z  =  W  c  +  v  (38) 

where  z,  c,  and  v  are  vectors  and  W  is  a  matrix  (as  are  all  boldface  quantities  in  what 
follows). 

We  can  define  the  correlation  matrix  of  the  received  signal 


R  =  <z  zH>  =  W  <c  cH>  WH  =  W  C  WH  +  o2  I  (39) 

where  brackets  denote  the  statistical  ensemble  average.  This  correlation  matrix  has  a 
special  structure,  the  exploitation  of  which  is  an  integral  part  oi  the  MUSIC  algorithm. 
First,  because  the  steering  vectors  are  deterministic  quantities,  the  brackets  pass  through 
them  to  the  (stochastic)  signal  amplitudes  on  the  inside,  thus  defining  the  power  matrix  C. 
If  the  signal  amplitudes  are  uncorrelated  among  themselves,  then  C  has  full  rank  K,  and  as 
we  will  see,  this  is  a  necessary  condition  for  MUSIC  to  work.  Since  the  noise  is  assumed 
uncorrelated  from  element  to  element,  the  brackets  give  an  N  by  N  (K<N)  diagonal  matrix 
for  the  correlation  of  the  noise  with  itself,  and  since  the  noise  is  assumed  stationary  over 
the  array  (as  well  as  over  the  interval  of  observation)  this  diagonal  matrix  is  simply  the 
identity  matrix  multiplied  by  a  scalar,  the  noise  power  Since  the  noise  is  uncorrelated 
with  the  signal,  there  are  no  cross  terms  between  noise  and  signal. 

The  correlation  matrix  R  satisfies  an  eigenvalue  equation 

R  =  j  =  0 . N-l  .  (40) 

These  eigenvalues  have  an  important  property,  namely,  that  N-K  of  them  are  degenerate 
and  equal  to  the  noise  power 


=  C'  j  =  K  +  l . N-l. 


.3 


(41) 
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Furthermore,  this  degenerate  eigenvalue  is  the  minimum  eigenvalue',  all  others  are  greater. 
The  argument  is  based  on  the  rank  structure  of  the  correlation  matrix  R.  The  signal  part  of 
this  matrix  (WCWh)  has  rank  K  (due  to  a  theorem  of  linear  algebra  that  says  that  the  rank 
of  a  product  of  three  matrices  is  controlled  by  the  one  in  the  middle),  and  so  this  term  can 
only  contribute  towards  K  of  the  eigenvalues.  The  noise  part  of  the  correlation  matrix, 
however,  has  full  rank  N,  and  so  contributes  to  all  eigenvalues.  (Being  <j~  times  the  unit 
matrix,  this  contribution  is  just  a2-)  Furthermore,  the  signal  part  of  the  correlation  matrix 
contributes  only  positive  amounts  to  the  eigenvalues  (in  addition  to  the  G2  from  the  norse), 
because  it  is  Hermitian  and  positive  definite.  (A  Hermitian  matrix  has  only  real 
eigenvalues,  and  a  positive  definite  matrix  has  positive  eigenvalues.)  So  the  noise 
eigenvalues  are  the  minimum  eigenvalues. 

With  this  assumption  on  the  eigenvalues,  we  can  write 


R^  =  WCW'h^  +  g2  =  o2  (j  =  K+1, .  .  .  ,  N-l) 

=>  \V  C  WH  =  0  (j=  K+1,...,  N-l) 

=>  WH  =  0  (j  =  K+1,...,  N-l) 

=>  WH  E  j  =  0  Ej  =  £j  O'  =  K+1,...  ,  N-l). 

The  next  to  last  step  is  true  only  if  the  matrix  C  has  full  rank  (so  that  it  is  nonsingular  and 
can  be  inverted),  and  since  this  is  the  key  to  the  MUSIC  method,  we  see  why  it  is 
important  for  the  signals  to  be  uncorrelated. 

We  can  therefore  decompose  the  eigenspace  into  a  K-dimensional  signal  subspace  spanned 
by  the  eigenvectors  associated  with  the  K  eigenvalues  that  are  greater  than  a2,  and  an  N-K 
dimensional  noise  subspace  of  eigenvectors  Ej  that  are  orthogonal  to  the  steering  vector  (as 
indicated  above). 

This  orthogonality  occurs  only  when  the  steering  vector  is  in  fact  steering  in  the 
direction  of  the  signal,  as  assumed  in  the  construction  above.  Therefore,  we  can  calculate 
or  estimate  the  steering  vector  for  a  whole  grid  of  directions  (this  is  sometimes  called  the 
array  manifold),  solve  the  eigenvalue  problem  of  the  correlation  matrix,  and  search  over  the 
array  manifold  for  the  directions  that  are  orthogonal  to  the  noise  eigenvectors.  To  do  this, 
we  must  first  estimate  the  correlation  matrix  from  the  data  samples.  Averaging  over 
snapshots  will  asymptotically  lead  to  the  true  correlation  matrix  (by  ergodicity),  but  of 
course  we  cannot  average  for  an  interval  longer  than  that  for  which  the  target  and  noise 
characteristics  remain  stationary.  The  averaged  data  correlation  matrix  has  the  form 

R  -  ±  I  Z(t,>  z»(t,) .  (43) 

M  j .  i 

A  MUSIC  pseudospectrum  can  be  generated  by  the  following  expression 


S  (0)  =  - — 1 - 

WH(0)  EE11  W(0) 

E  -  (Ek+i  •  •  •  En  i  I  • 


(44) 


Note  that  sharp  peaks  occur  when  the  denominator  is  zero,  i.e  ,  when  the  orthogonality 
condition  is  reached.  The  denominator  is  actually  (the  square  of)  the  average  over  all  noise 
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eigenvectors  of  the  equation  WH  Ej  =  0.  This  reduces  the  effect  of  spurious  zeros  that  do 
not  represent  a  true  signal  but  that  may  arise  when  the  dimension  of  the  noise  space  is 
greater  than  the  number  of  signals.  The  true  zeros  persist  over  the  averaging  because  they 
satisfy  the  orthogonality  condition,  but  the  accidental  zeros  do  not. 

3.  MUSIC  RESULTS 

The  MUSIC  results  include  all  figures  beginning  with  the  letters  "A"  and  "B."  They 
are  organized  as  follows. 

All  cases  have  1 1  elements  with  XJ2  spacing.  All  cases  have  two  equal  amplitude  plane 
waves  coming  in:  one  at  Q\  =  45°  and  the  other  at  02  =  45°  to  50°  at  1/2°  spacing.  (All 
figures  show  the  true  angles  of  arrival  as  straight  lines  on  the  lower  subplots.)  Figures 
beginning  with  the  letter  "A"  all  have  perfect  decorrelation  between  the  signals  (a  =  1.0), 
and  figures  beginning  w'ith  the  letter  ”B”  consider  cases  at  or  near  perfect  correlation  (a  = 
0.5  for  perfect  correlation). 

Note:  The  signals  are  generated  as  follows  (for  X/2  spacing): 

Si(JSNAP,  IELEM)  =  ej(JSNAP'I>i7  K  sin  01 

S2(JSNAP,  IELEM)  =  e-1  a(JSNAP-l)3.  ^(IELEM-l)  n  sin  02 

S total  =  Si  +  S2  +  v  (JSNAP, IELEM) 

where  Sj  is  the  first  plane  wave  for  snapshot  JSNAP  at  element  IELEM,  S2  is  the  second 
plane  wave,  and  v  is  the  white  Gaussian  uncorrelated  noise*  .  In  this  case,  for  a  =  1.0, 
the  signals  ("in  space,"  see  Ref.  5)  are  uncorrelated  due  to  Fourier  orthogonality. 
(Statistical  orthogonality,  i.e.,  decorrelation,  is  a  much  more  general  stochastic  concept  that 
includes  Fourier  orthogonality  as  a  particular  case,  at  least  for  the  ergodic  case  where 
ensemble  averages  are  equal  to  integrals  over  time.) 

3.1  Explanation  of  MUSIC  Figures 

The  "A"  group  is  divided  into  four  series,  which  we  will  refer  to  as  A 1  (Figs.  A  1.1  to 
Al.ll),  A2  (Figs.  A2.1  to  A2.8),  A3  (Figs.  A3.1  to  A3. 8),  and  A4  (Figs.  A4.1  to  A4.8). 

•  Series  A1  considers  the  case  of  64  snapshots  for  various  SNRs.  The  SNRs 
considered  are:  5,  7,  8,  9,  10,  15,  20,  30,  40,  50,  and  99  dB. 

•  Series  A2  considers  the  case  of  10  dB  SNR  for  various  numbers  of  snapshots.  The 
numbers  of  snapshots  considered  are:  1,  2,  5,  10,  20,  30,  64,  and  100  snapshots. 

•  Series  A3  considers  the  case  of  20  dB  SNR  for  various  numbers  of  snapshots.  The 
numbers  of  snapshots  considered  are:  1,  2,  5,  10,  20,  30,  64,  and  100  snapshots. 


(45) 

(46) 

(47) 


The  noise  is  generated  according  to  a  standard  algorithm  for  converting  uniformly  distributed  random  variables 
(as  provided  by  the  random  number  generator  of  a  computer)  to  while  Gaussian  random  variables.  A  and  B  are 
random  variables  uniformly  distributed  in  [0.1).  The  variance  s^  is  derived  (for  unity  amplitude  signals)  as 


a2  = 


1 


210" 


With  R  =  ,|2<73  In 


,  \-A 


and  0  —  2/rfi  a  sample  of  complex  white  Gaussian  noise  v  is  given  by  Re  ‘H. 
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•  Series  A4  considers  the  case  of  30  dB  SNR  for  various  numbers  of  snapshots.  The 
numbers  of  snapshots  considered  are:  1,  2,  5,  10,  20,  30,  64,  and  100  snapshots. 

The  "B"  group  is  divided  into  two  series,  which  we  will  refer  to  as  B1  (Figs.  Bl.l  to 
B1.7)  and  B2  (Figs.  B2.1  to  B2.7). 

•  Series  B1  considers  the  case  of  64  snapshots  and  10  dB  SNR  for  various  a's.  The 
a's  considered  are:  0.50  (p  =  1.00),  0.55  (p  =  0.94),  0.56  (p  =  0.91),  0.57  (p  = 
0.88),  0.58  (p  =  0.84),  0.59  (p  =  0.80),  and  0.60  (p  =  0.76). 

•  Series  B2  considers  the  case  of  64  snapshots  and  20  dB  SNR  for  various  a’s.  The 
a's  considered  are:  0.50  (p  =  1.00),  0.51  (p  =  1.00),  0.52  (p  =  0.99),  0.53  (p  - 
0.98),  0.54  (p  =  0.96),  0.55  (p  =  0.94),  and  0.56  (p  =  0.91). 

All  MUSIC  figures  have  the  same  format:  an  upper  plot  with  a  MUSIC 
pseudospectrum,  and  a  lower  plot  representing  the  estimates  of  the  angles  of  arrival.  Note 
that  the  upper  plot  also  contains  the  tops  of  the  broad  peaks  resulting  from  an  unwindowed 
FFT  of  the  datapoints  in  the  first  of  the  snapshots  from  the  set  of  snapshots  that  generated 
the  corresponding  MUSIC  spectrum.  This  is  to  indicate  the  dramatic  improvement  in 
resolution  permitted  by  the  high-resolution  algorithms.  The  estimates  of  the  angles  of 
arrival  in  the  lower  plot  were  found  simply  by  searching  for  local  maxima  in  the  MUSIC 
pseudospectrum.  The  resolution  of  the  MUSIC  spectrum  is  256  points  over  the  angular 
range  in  the  axis  of  the  spectral  plot.  The  true  angles  of  arrival  are  indicated  by  the  straight 
lines  in  the  lower  plot. 

3.1.1  Series  A1  (a  =  1.00;  64  snaps;  vancus  SNRs) 

The  general  trend  is  that  the  signals  are  resolved  at  an  increasingly  lower  angular 
separation  with  increasingly  high  SNR,  until  at  50  dB  SNR  and  above,  the  estimates  are 
essentially  perfect.  There  is  bias  in  the  estimates  near  the  angular  threshold  of  resolution, 
as  one  might  expect  (one  apparent  peak  gradually  splits  into  two),  but  the  estimates  rapidly 
converge  to  the  true  values  above  this  threshold.  Interestingly,  the  lower  angle  (i.e.,  45°) 
quickly  converges  to  its  exact  value,  while  the  upper  angle  keeps  a  slight  bias  over  the 
angular  range  indicated  until  the  SNR  reaches  40  dB.  This  asymmetry  has  nothing  to  do 
with  the  fact  that  the  lower  angle  is  held  constant,  because  the  estimates  are  made  separately 
for  the  different  angular  separations.  The  resolution  thresholds  for  the  various  SNRs  are: 


-  8  dB  SNR: 

P*  =  0.27  (4.6°)  (of  angular  separation) 

-  9  dB  SNR: 

P  =  0.24  (4.0°) 

-  10  dB  SNR 

P  =  0.24  (4.0°) 

-  15  dB  SNR 

P  =0.18  (3.0°) 

-  20  dB  SNR 

P  =  0.12  (2.0°) 

-  30  dB  SNR 

P  =0.06  (1.0°) 

-  40  dB  SNR 

P  =  0.03  (0.5°) 

-  50  dB  SNR 

p  =  0.00  (0.0°) 

-  99  dB  SNR 

p  =  0.00  (0.0°). 

3.1 .2  Series  A2  (a  = 

1.00;  10  dB  SNR;  various  numbers  of  snapshots) 

Resolution  is  not  attained  (over  the  angular  range  of  these  plots)  until  20  snapshots. 
With  more  snapshots,  the  resolution  threshold  tends  to  decrease,  although  at  64  snapshots 
it  increased  over  the  previous  case,  showing  that  more  snapshots  does  not  necessarily 
result  in  an  improvement.  At  30  snapshots,  the  lower  angle  rapidly  converges  to  its  true 

*(3  is  a  measure  of  angular  separation  of  targets  that  takes  into  account  the  reduced  effective  aperture  of  the  array  as 
wc  scan  away  from  boresight:  P  =  5  |sin02  -  sinflj ). 
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value,  but  the  upper  angle  maintains  a  slight  bias,  and  this  pattern  is  repeated  even  at  100 
snapshots.  This  parallels  the  trend  observed  in  Series  Al,  in  which  the  upper  angle  keeps  a 
bias  (until  high  SNR)  while  the  lower  angle  does  not.  The  resolution  thresholds  for  the 
various  numbers  of  snapshots  are: 

-  20  snapshots:  p  =  0.24  (4.0°) 

-  30  snapshots:  (3=0.21  (3.5°) 

-  64  snapshots:  P  =  0.24  (4.0°) 

-  100  snapshots:  P  =0.18  (3.0°). 

3.1 .3  Series  A3  (a  =  1 .00;  20  dB  SNR;  various  numbers  of  snapshots) 

Resolution  is  not  attained  (over  the  angular  range  of  these  plots)  until  10  snapshots. 
With  more  snapshots,  the  resolution  threshold  tends  to  decrease,  at  first  quickly,  then  much 
more  slowly.  Again,  the  upper  angle  tends  to  have  slightly  more  bias  over  all  cases  than 
the  lower  angle,  although  we  can  say  that  the  estimates  are  very  good  above  the  threshold 
for  20  or  more  snapshots.  The  resolution  thresholds  for  the  various  numbers  of  snapshots 
are: 

-  10  snapshots:  P  =  0.27  (4.5°) 

-  20  snapshots:  P  =  0.15  (2.5°) 

-  30  snapshots:  P  =  0.12  (2.0°) 

-  64  snapshots:  P  =  0.12  (2.0°) 

-  100  snapshots:  P  =  0.09  (1,5°). 

3.1.4  Series  A4  (a  =  1.00;  30  dB  SNR;  various  numbers  of  snapshots) 

Resolution  is  attained  (over  the  angular  range  of  these  plots)  at  just  5  snapshots.  With 
more  snapshots,  the  resolution  threshold  tends  to  decrease,  at  first  quickly,  then  more 
slowly.  The  upper  angle  no  longer  tends  to  have  slightly  more  bias  than  the  lower  angle. 
The  estimates  are  good  for  10  snapshots  and  very  good  for  20  or  more  snapshots,  with 
rapid  convergence  above  threshold  to  the  true  values.  The  resolution  thresholds  for  the 
various  numbers  of  snapshots  are: 

-  5  snapshots:  P  =  0.24  (4.0°) 

-  10  snapshots:  P  =  0. 15  (2.5°) 

-  20  snapshots:  P  =  0.06  (1.0°) 

-  30  snapshots:  P  =  0.06  (1.0°) 

-  64  snapshots:  P  =  0.06  (1.0°) 

-  100  snapshots:  p  =  0.03  (0.5°). 

3.1 .5  Series  B1  (64  snapshots;  10  dB  SNR;  various  as) 

Resolution  is  attained  (over  the  angular  range  of  these  plots)  at  a  =  0.57  (p  =  0.88). 
The  resolution  threshold  decreases  somewhat  with  higher  a  (although  we  do  not  consider 
a's  higher  than  0.60).  The  resolution  thresholds  for  the  various  a's  are: 

-  a  =  0.57  (p  =  0.88):  p  =  0.27  (4.5°) 

-  a  =  0.58  (p  =  0.84):  p  =  0.24  (4.0°) 

-  a  =  0.59  (p  =  0.80):  p  =  0.24  (4.0°) 

-  a  =  0.60  (p  =  0.76):  p  =  0.24  (4.0°). 
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3.1 .6  Series  B2  (64  snapshots;  20  dB  SNR;  various  a's) 

Resolution  is  attained  (over  the  angular  range  of  these  plots)  at  a  =  0.53  (p  =  0.98). 
The  resolution  threshold  decreases  somewhat  with  higher  a  (although  we  do  not  consider 
a's  higher  than  0.56).  The  resolution  thresholds  for  the  various  a's  are: 


-  a  =  0.53  (p  =  0.98): 

(3  =  0.24 

(4.0°) 

-  a  =  0.54  (p  =  0.96): 

P  =  0.21 

(3.5°) 

-  a  =  0.55  (p  =  0.94): 

P  =  0.21 

(3.5°) 

-  a  =  0.56  (p  =0.91): 

P  =0.18 

(3.0°). 

3.2  MUSIC:  Conclusions 

3.2.1  Perfectly  uncorrelated  case 

The  resolution  properties  of  the  MUSIC  algorithm  are  excellent  in  the  perfectly 
uncorrelated  case  (a  =  1.0).  The  resolution  threshold  decreases  with  both  increased  SNR 
and  increased  number  of  snapshots,  with  SNR  being  the  more  important  factor.  For  64 
snapshots,  a  resolution  of  (3  =  0.12  (2.0°)  is  achieved  at  20  dB  SNR  and  of  (3  =  0.06 
(1.0°)  at  30  dB  SNR.  Above  50  dB  SNR,  the  estimates  are  essentially  perfect.  For  10  dB 
SNR,  resolution  is  attained  at  (3  =  0.24  (4.0°)  with  20  snapshots.  For  20  dB  SNR, 
resolution  is  attained  with  10  snapshots,  but  is  much  better  with  20  or  more  snapshots 
(1.5°  to  2.5°).  For  30  dB  SNR,  resolution  is  attained  with  5  snapshots,  but  is  much  better 
with  10  or  more  snapshots  (2.5°  for  10  snapshots;  0.5°  to  1.0°  for  20  or  more  snapshots). 
Above  the  resolution  threshold,  the  estimates  rapidly  converge  to  the  true  values,  especially 
with  20  or  more  snapshots  at  SNRs  of  20  dB  or  higher.  At  lower  SNR  and  fewer 
snapshots,  there  is  a  tendency  towards  more  bias  in  the  upper  angle  than  in  the  lower  one 
(for  the  angular  cases  considered  here). 

3.2.2  Nearly  correlated  cases 

The  MUSIC  algorithm  can  achieve  resolution  at  surprisingly  high  correlation, 
considering  the  theoretical  (and  actual)  breakdown  at  perfect  correlation.  For  10  dB  SNR 
and  64  snapshots,  resolution  is  attained  at  a  =  0.57  (p  =  0.88),  although  the  angular 
separation  (4.5°)  is  probably  too  large  to  be  considered  "high  resolution."  (This  separation 
does  not  decrease  as  a  approaches  0.60  (p  =  0.76),  the  maximum  considered  here.)  For 
20  dB  SNR  and  64  snapshots,  resolution  is  attained  at  a  =  0.53  (p  =  0.98),  with  an 
angular  separation  of  4.0°.  This  separation  decreases  to  3.0°  as  a  approaches  0.56  (p  = 
0.91),  the  maximum  considered  here. 

4.  PRONY  RESULTS 

The  Prony  results  include  all  Figures  beginning  with  the  letter  "C."  They  are  organized 
as  follows. 

All  cases  have  1 1  elements  with  X/2  spacing.  All  cases  have  two  equal  amplitude  plane 
waves  coming  in:  one  at  0j  =  45°  and  the  other  at  02  =  45°  to  60°  at  1/2°  spacing.  (All 
Figures  show  the  true  angles  of  arrival  as  straight  lines.)  Most  importantly,  all  cases 
consider  the  perfect  correlation  case  (a  =  0.5),  since,  as  we  have  seen,  such  eigenvector 
techniques  as  MUSIC  work  very  well  when  just  a  little  decorrelation  is  present.  It  is  in  the 
case  of  perfect  correlation  (which  may  occur,  for  example,  with  multipath)  that  makes  the 
Prony  technique  particularly  interesting.  We  consider  both  the  "Prony  space"  and  the 
"Prony  space-time"  variations  of  the  algorithm,  as  described  previously  in  this  report. 
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4.1  Explanation  of  Prony  Figures 

There  are  four  series,  which  we  will  refer  to  as:  Cl  (Figs.  Cl.l  to  Cl. 9),  C2  (Figs. 
C2.1  to  C2.9),  C3  (Figs.  C3.1  to  C3.9),  and  C4  (Figs.  C4.1  to  C4.9).  Within  a  series, 
the  parameter  that  varies  is  the  number  of  snapshots.  We  will  consider  the  following  cases: 
1,  2,  3,  5,  10,  20,  30,  64,  and  100  snapshots.  From  series  to  series,  the  parameter  that 
varies  is  the  SNR  (defined  previously  in  this  report).  Series  Cl  considers  10  dB  SNR;  C2 
considers  20  dB  SNR;  C3  considers  30  dB  SNR;  and  C4  considers  40  dB  SNR.  in  this 
way,  we  have  a  Cartesian  product  of  SNR  vs  number  of  snapshots,  with  all  other 
parameters  maintained  constant.  Our  objective  is  to  resolve  the  signals  with  the  fewest 
snapshots  for  a  given  SNR. 

4.1.1  Series  Cl  (10  dB  SNR) 

In  all  cases,  the  Prony  space  algorithm  shows  a  threshold  at  about  14°  of  angular 
separation  between  the  plane  waves  (0i  =  45°  and  02  =  59°).  This  corresponds  to  p  = 
0.75.  The  threshold  suddenly  becomes  sharp  at  3  snapshots  (Fig.  Cl. 3),  although  we  can 
say  that  the  signals  are  resolved  (with  bias)  at  1  and  2  snapshots.  With  1  snapshot  (Fig. 
Cl.l)  the  bias  at  P  =  0.75  is  a  few  degrees,  and  with  2  snapshots  (Fig.  Cl. 2)  no  bias 
exists  for  Prony  space  at  the  upper  angle  (02  =  59°)  and  a  couple  of  degrees  of  bias  exist  at 
the  lower  angle  (0i  =  45°).  At  3  snapshots  (Fig.  Cl. 3),  very  little  bias  (at  P  =  0.75)  is 
present  on  the  lower  angle,  but  there  is  a  bias  of  about  5  degrees  at  the  upper  angle,  which 
decreases  appreciably  with  an  increase  in  the  number  of  snapshots. 

The  Prony  space-time  algorithm  works  slightly  better,  with  about  the  same  threshold  of 
resolution  (p  =  0.75),  but  without  the  bias  on  the  upper  angle  at  the  threshold  itself.  In 
fact,  one  could  say  that  the  Prony  space-time  algorithm  generally  resolves  the  signals  at  a 
much  lower  threshold  (P  @  0.4  or  0.5),  but  with  erratic  and  often  large  biases  below  the 
threshold  of  the  Prony  space  algorithm. 

4.1.2  Series  C2  (20  dB  SNR) 

In  all  cases,  the  Prony  space  algorithm  shows  a  threshold  at  about  8°  of  angular 
separation  between  the  plane  waves  (0i  =  45°  and  02  =  53°).  This  corresponds  to  P  = 
0.46.  The  threshold  suddenly  becomes  sharp  at  3  snapshots  (Fig.  C2.3),  although  we  can 
say  that  the  signals  are  resolved  (with  bias)  at  1  and  2  snapshots.  With  1  snapshot  (Fig. 
C2.1)  the  bias  at  p  =  0.46  is  a  few  degrees,  and  with  2  snapshots  (Fig.  C2.2)  there  is  no 
bias  for  Prony  space  at  the  upper  angle  (02  =  53°)  and  a  couple  of  degrees  of  bias  at  the 
lower  angle  (0i  =  45°).  At  3  snapshots  (Fig.  C2.3),  there  is  no  bias  (at  P  =  0.46)  on  the 
lower  angle,  but  a  bias  of  about  7  degrees  at  the  upper  angle,  which  decreases  appreciably 
with  an  increase  in  the  number  of  snapshots.  Above  the  threshold,  the  Prony  space 
estimates  rapidly  converge  to  the  true  values. 

The  Prony  space-time  algorithm  works  slightly  better,  with  about  the  same  threshold  of 
resolution  (p  =  0.46),  but  without  such  a  large  bias  on  the  upper  angle  at  the  threshold 
itself.  In  fact,  one  could  say  that  the  Prony  space-time  algorithm  generally  resolves  the 
signals  at  a  much  lower  threshold  (p  @  0.35  or  even  0.25  in  some  cases),  but  with  erratic 
and  often  large  biases  below  the  threshold  of  the  Prony  space  algorithm.  (The  pattern  of 
these  erratic  biases  bears  some  similarity  to  those  of  the  previous  series.) 

4.13  Series  C 3  (30  dB  SNR) 

In  all  cases,  the  Prony  space  algorithm  shows  a  threshold  at  about  4.5°  of  angular 
separation  between  the  plane  waves  (0i  =  45°  and  02  =  49.5°).  This  corresponds  to  p  = 
0.27.  The  threshold  suddenly  becomes  sharp  at  3  snapshots  (Fig.  C3.3),  although  we  can 
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say  that  the  signals  are  resolved  (with  bias)  at  1  and  2  snapshots.  With  1  snapshot  (Fig. 
C3.1),  the  bias  at  P  =  0.27  is  a  few  degrees  or  less,  and  with  2  snapshots  (Fig.  C3.2), 
there  is  negligible  bias  for  Prony  space  at  the  upper  angle  (02  =  49.5°)  and  only  a  degree  of 
bias  at  the  lower  angle  (0i  =  45°).  At  3  snapshots  (Fig.  C3.3),  there  is  no  bias  (at  (3  = 
0.27)  on  the  lower  angle,  but  a  bias  of  about  5  degrees  at  the  upper  angle,  which  decreases 
appreciably  with  an  increase  in  the  number  of  snapshots.  Above  the  threshold,  the  Prony 
space  estimates  rapidly  converge  to  the  true  values. 

The  Prony  space-time  algorithm  works  slightly  better,  with  about  the  same  threshold  of 
resolution  (f3  =  0.27),  but  without  such  a  large  bias  on  the  upper  angle  at  the  threshold 
itself.  In  fact,  one  could  say  that  the  Prony  space-time  algorithm  generally  resolves  the 
signals  at  a  much  lower  threshold  (P  @  0.20  or  even  0.15  in  some  cases),  but  with  erratic 
and  often  large  biases  below  the  threshold  of  the  Prony  space  algorithm.  (The  pattern  of 
these  erratic  biases  bears  some  similarity  to  those  of  the  two  previous  series.) 

4.1 .4  Series  C4  (40  dB  SNR) 

In  all  cases,  the  Prony  space  algorithm  shows  a  threshold  at  about  2.5°  of  angular 
separation  between  the  plane  waves  (0j  =  45°  and  02  =  47.5°).  This  corresponds  to  P  = 
0.15.  The  threshold  suddenly  becomes  sharp  at  3  snapshots  (Fig.  C4.3),  although  we  can 
say  that  the  signals  are  resolved  (with  bias)  at  1  and  2  snapshots.  With  1  snapshot  (Fig. 
C4.1),  the  bias  at  p  =  0.15  is  a  few  degrees  or  less,  and  with  2  snapshots  (Fig.  C3.2), 
there  is  negligible  bias  for  Prony  space  at  the  upper  angle  (02  =  47.5°)  and  only  a  degree  of 
bias  at  the  lower  angle  (0i  =  45°).  At  3  snapshots  (Fig.  C4.3),  there  is  no  bias  (at  P  = 
0.15)  on  the  lower  angle,  but  a  bias  of  about  5  degrees  at  the  upper  angle,  which  decreases 
appreciably  with  an  increase  in  the  number  of  snapshots.  Above  the  threshold,  the  Prony 
space  estimates  rapidly  converge  to  the  true  values. 

The  Prony  space-time  algorithm  works  slightly  better,  with  about  the  same  threshold  of 
resolution  (P  =  0.15),  but  without  such  a  large  bias  on  the  upper  angle  at  the  threshold 
itself.  In  fact,  one  could  say  that  the  Prony  space-time  algorithm  generally  resolves  the 
signals  at  a  much  lower  threshold  (P  @  0.10  or  even  0.09),  but  with  erratic  and  often  large 
biases  below  the  threshold  of  the  Prony  space  algorithm. 

4,2  Prony:  Conclusions 

The  different  series  of  cases  (different  SNRs)  just  considered  obviously  exhibit  a 
similar  pattern.  The  Prony  space  algorithm  show  a  well-defined  threshold  for  good  (i.e., 
unbiased)  resolution  that  decreases  with  SNR:  P  =  0.75  (14°  of  angular  separation)  at  10 
dB  SNR;  P  =  0.46  (8°  of  angular  separation)  at  20  dB  SNR;  p  =  0.27  (4.5°  of  angular 
separation)  at  30  dB  SNR;  and  P  =  0.15  (2.5°  of  angular  separation)  at  40  dB  SNR.  This 
threshold  always  becomes  sharp  at  just  3  snapshots,  although  in  all  cases  we  can  say  that 
we  have  resolved  the  signals  (with  bias)  at  just  1  snapshot.  At  3  snapshots,  there  is  always 
a  bias  of  about  5°  (at  threshold)  at  the  upper  angle,  which  decreases  appreciably  with  an 
increase  in  the  number  of  snapshots,  and  there  is  little  bias  (at  threshold)  at  the  lower  angle. 
For  20  dB  or  more  of  SNR,  the  estimates  converge  rapidly  above  the  threshold  to  the  true 
values,  and  even  at  10  dB  SNR,  the  biases  are  not  large  above  the  threshold.  The  Prony 
space-time  algorithm  works  at  least  as  well  as  the  Prony  space  algorithm,  and  if  one  is 
willing  to  tolerate  erratic  and  large  biases  near  threshold,  one  can  say  that  the  Prony  space- 
time  algorithm  shows  an  improvement  in  the  P-threshold  of  0.1  or  even  0.2  over  the  Pronv 
space  algorithm  (with  the  larger  improvement  at  lower  SNR). 

Actually,  it  appears  that  with  3  snapshots,  resolution  (with  bias)  occurs  in  Prony  space- 
time  down  to  p=0  at  all  SNR’s  (10,  20,  30,  and  40  dB).  Although  the  bias  is  not  small, 
this  result  may  still  be  significant  in  that  at  least  it  tells  us  that  another  target  is  there,  in 
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approximately  a  certain  direction.  One  must  still  ask,  though,  whether  this  ’resolution'  is 
spurious.  Why,  for  instance,  does  it  disappear  with  more  snapshot  averaging?  Even 
though  this  phenomenon  is  curious,  it  seems  difficult  to  believe  that  it  is  completely 
accidental. 

5.  CONCLUSIONS 

Detailed  conclusions  for  the  MUSIC  and  Prony  algorithms  were  provided  in  Sections  3 
and  4,  respectively. 

MUSIC  gives  excellent  resolution  in  all  but  the  most  correlated  cases.  Following  the 
advice  of  the  reviewer,  we  present  in  Appendix  B  some  results  from  a  variation  of  MUSIC 
called  spatially  averaged  MUSIC  that  performs  as  well  as  Prony  in  the  case  of  perfectly 
correlated  signals.  The  resolution  threshold  decreases  with  both  increased  SNR  and 
increased  number  of  snapshots,  with  SNR  being  the  more  important  factor.  In  the 
perfectly  uncorrelated  case,  and  for  64  snapshots,  a  resolution  of  p  =  0.12  (2.0°)  is 
achieved  at  20  dB  SNR  and  of  (3  =  0.06  (1.0°)  at  30  dB  SNR.  Above  50  dB  SNR,  the 
estimates  are  essentially  perfect.  For  the  nearly  correlated  case  of  a  =  0.53  (perfect 
correlation:  a  =  0.50),  an  angular  separation  of  4.0°  is  attained  for  20  dB  SNR  and  64 
snapshots. 

For  the  Prony  space  algorithm  in  the  perfectly  correlated  case,  there  is  a  well-defined 
threshold  for  good  (i.e.,  unbiased)  resolution  that  decreases  with  SNR:  14°  of  angular 
separation  at  10  dB  SNR;  8°  of  angular  separation  at  20  dB  SNR;  4.5°  of  angular 
separation  at  30  dB  SNR;  and  2.5°  of  angular  separation  at  40  dB  SNR.  This  threshold 
always  becomes  sharp  at  just  3  snapshots,  although  in  all  cases  we  can  say  that  we  have 
resolved  the  signals  (with  bias)  at  just  1  snapshot.  Prony  space-time  seems  to  work  even 
better.  In  some  cases,  the  ’’resolution"  is  so  good  that  one  wonders  if  it  is  real. 

In  summary,  the  well-known  excellent  performance  of  non-spatially-averaged  MUSIC 
in  all  but  the  most  correlated  cases  has  been  verified,  and  it  has  been  demonstrated  that  the 
Prony  technique  is  promising  in  the  perfectly  correlated  case.  Also,  the  Prony  technique 
seems  to  require  few  snapshots  to  resolve  close  targets. 
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Rngular  Separation  (deg)  Angular  Separation  (deg) 

Fig.  Al.l  -  MUSIC  algorithm;  perfect  decorrelation;  64  snapshots;  SdBSNR  Fig.  A  1.2-  MUSIC  algorithm;  perfect  decorrelation,  64  snapshots,  7dBSNR 


NRL  REPORT  9397 


19 


Angular  Separation  (deg)  Rngulor  Separation  (deg) 

Fig.  A1.5  -  MUSIC  algorithm;  perfect  decorrelation;  64  snapshots;  10  dB  SNR  Fig.  A1.6  -  MUSIC  algorithm;  perfect  decorrelation;  64  snapshots;  13  dB  SNR 


Angular  Separation  (deg)  Rngulor  Separation  (deg) 

Fig.  A4.1-  MUSIC  algorithm;  perfect  decorrelation;  1  snapshot;  30dBSNR  Fig.  A4.2-  MUSIC  algorithm;  perfect  decorrelation;  2  snapshots;  30dBSNR 
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Fig.  B  1.5-  MUSIC  algorithm;  high  correlation  (a=0.58);  64  snapshots;  10  dB  SNR  Fig.  B  1.6-  MUSIC  algorithm;  high  correlation  (a=0.59);  64  snapshots;  10  dB  SNR 


fingular  Separation  (deg)  Rngulor  Separation  (deg) 

Fig.  B2.1-  MUSIC  algorithm;  high  correlation  (a=0.50);  64  snapshots;  20dBSNR  Fig.  B2.2-  MUSIC  algorithm;  high  correlation  (a=0.51);  64  snapshots;  20dBSNR 
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flngulor  Separation  (deg)  Angular  Separation  (deg) 

Fig.  Cl. I  -  PRONY  algorithm;  perfect  correlation;  10  dB  SNR;  1  snapshot  Fig.  Cl. 2  -  PRONY  algorithm;  perfect  correlation;  10  dB  SNR;  2  snapshots 
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flng-jLar  Separation  (deg)  finqular  Separation  (deg) 

Fig.  Cl. 3  -  PRONY  algorithm;  perfect  correlation;  10  dB  SNR;  3  snapshots  Fig.  Cl. 4  -  PRONY  algorithm;  perfect  correlation;  10  dB  SNR;  5  snapshots 
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Angular  Separation  (deg)  Angular  Separation  (deg) 

Fig.  Cl  .5  -  PRONY  algorithm;  perfect  correlation;  10  dB  SNR;  10  snapshots  Fig.  Cl  .6  -  PRONY  algorithm;  perfect  correlation;  10  dB  SNR;  20  snapshots 
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Angular  Separation  (deg) 

Fig.  Cl. 9  -  PRONY  algorithm;  perfect  correlation;  10  dB  SNR;  100  snapshots 
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RngyLor  SeporotLon  (deg)  Angular  Separation  (deg) 

I-ig.  C2.1  -  PRONY  algorithm;  perfect  correlation;  20  dB  SNR;  1  snapshot  Fig.  C2.2  -  PRONY  algorithm;  perfect  correlation;  20  dB  SNR;  2  snapshots 
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flngulor  Separation  (deg)  Angular  Separation  ( deq ) 

Fig.  C2.5  -  PRONY  algorithm;  perfect  correlation;  20  dB  SNR;  10  snapshots  Fig.  C2.6  -  PRONY  algorithm;  perfect  correlation;  20  dB  SNR;  20  snapshots 
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Angular  Separation  (deg)  flngulor  Separation  (deg) 

Fig.  C2.7  -  PRONY  algorithm;  perfect  correlation;  20  dB  SNR;  30  snapshots  Fig.  C2.8  -  PRONY  algorithm;  perfect  correlation;  20  dB  SNR;  64  snapshots 
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Angular  Separation  (deg)  Angular  Separation  (deg) 

Rfr  C3. 1  -  PRONY  algorithm;  perfect  correlation;  30  dB  SNR;  1  snapshot  Fig.  C3.2  -  PRONY  algorithm;  perfect  correlation;  30  dB  SNR;  2  snapshots 
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Fig.  C3.3  -  PRONY  algorithm;  perfect  correlation;  30  dB  SNR;  3  snapshots  Fig.  C3.4  -  PRONY  algorithm;  perfect  correlation;  30  dB  SNR;  5  snapshots 
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Angular  Separation  (deg) 

Fig.  C3.9  -  PRONY  algorithm;  perfect  correlation;  30  dB  SNR;  100  snapshots 
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Angular  Separation  (deg)  Angular  Separation  (deg) 

Fig.  C4.3  -  PRONY  algorithm;  perfect  correlation;  40  dB  SNR;  3  snapshots  Fig.  C4.4  -  PRONY  algorithm;  perfect  correlation;  40  dB  SNR;  5  snapshots 
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Angular  Separation  (deg)  Angular  Separation  (deg) 

Fig.  C4.5  -  PRONY  algorithm;  perfect  correlation;  40  dB  SNR;  10  snapshots  Fig.  C4.6  -  PRONY  algorithm;  perfect  correlation;  40  dB  SNR;  20  snapshots 
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Angular  Separation  (deq) 

Fig.  C4.9  -  PRONY  algorithm;  perfect  correlation;  40  dB  SNR;  100  snapshots 


Appendix  A 

ESPRIT  ALGORITHM 


The  ESPRIT  algorithm  [Al]  is  an  attempt  to  achieve  performance  as  good  as  that  of 
MUSIC,  without  having  to  generate  the  time-  and  memory-consuming  array  manifold 
(steering  vectors).  ESPRIT  uses  two  or  more  arrays  that  are  translated  (or  rotated)  with 
respect  to  each  other,  and  then  takes  advantage  of  the  invariance  of  the  signal  subspace 
under  translation  to  solve  a  generalized  eigenvalue  equation  whose  eigenvalues  give  the 
directions  of  arrival.  (We  must  assume  stationarity  of  the  noise  between  the  arrays, 
although  we  do  not  need  to  assume  a  specific  structure  for  the  noise  component  of  the 
correlation  matrix,  e.g.,  diagonal,  as  we  did  with  MUSIC.)  In  the  following,  we  consider 
only  the  case  of  a  translated  array. 

The  signal  at  the  elements  of  the  shifted  array  is  related  to  that  at  the  elements  of  the 
unshifted  array  by  a  phase  change  [A2]: 

Ai  e'ikxx  — >  Ai  e"‘kx(x  +  A)  =  Ai  e_ikxx  e'ikxA  (Xm  =  m  d  ,  8  =  A  /  d) 

(Al) 

Al  e-kim  ->  Ai  e‘iki(m  +  5)  =  Aj  e^i™  c  lki5  (ki  =  kxd) . 

(This  is  the  projection  of  one  plane  wave  over  the  elements;  we  take  the  sum  of 
expressions  like  this  to  get  the  full  signal.  Ai  is  the  amplitude;  kx  is  the  projected 
wavevector;  x  is  the  position  along  the  array  in  real  space;  A  is  the  shift  of  the  second  array 
in  real  space;  m  is  the  number  of  the  element;  d  is  the  element  spacing;  ki  is  a 
dimensionless  wavevector  that  gives  the  signal  in  "digital"  space,  i.e.,  we  are  now 
analyzing  a  set  of  discrete  samples  with  normalized  frequencies  and  dimensionless  time 
interval  equal  to  1,  etc.) 

With  the  following  steering  vector 


k  =  2-^-sin(9)  S  =  [skl  sk2--- skL]  ,  (A2) 


we  can  write  for  the  total  signal  across  elements 

y  (n)  =  i  Ai(n)fiki  +  n(n)  .  (A3) 

i  =  1 


sk  = 


>ik 


>i2k 


e'Mkj 
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We  define  the  correlation  matrix  as  in  MUSIC 

R  =  E  [y(n)*  y(n)T]  =  £  sk.  Pjj  skj  +  o^,  I  (A4) 

i.j  =  1 

Pij  =  E  [  Ai(n)*  Aj(n)]  (l<i,j<L)  (A5) 

with  the  following  restrictions  on  the  noise  (although  the  left-hand  "white”  noise  condition 

is  not  necessary  in  the  general  case  if  the  noise  statistics  are  somehow  known) 

E[vm  (n)*  vk  (n)]  =  5mk  E[vm  (n)*  A,  (n)J  =  0  (A6) 

R  =  S  P  SH  +  I  .  (A7) 

The  steering  matrix  S  transforms  as 

e-‘ki5  0  0  0 

0  e'ik25  0  •  ■  •  0 

S  ->  S  D5  Dg  =  .  (A8) 


[  0  0  0  •••  e-ikLSj 

This  yields:  for  the  signal  at  the  unshifted  and  shifted  arrays, 

y(n)  =  S*  A(n)  +  v  (n)  (A9) 

y5(n)  =  S*  D*5  A(n)  +  v^n) ;  (A  10) 

for  the  unshifted  and  shifted  correlation  matrices, 

Ryy  =  E[y(n)*  y(n)T]  =  S  P  SH  +  I  (All) 

and 

Ry5y§  =  E[Xs(n)*  y§(n)T]  =  S  Dg  P  Dg  SH  +  I ;  (A  12) 

and  for  the  cross-correlation  matrix  between  the  unshifted  and  shifted  arrays, 

Ryyg  =  E  [y(n)*  y.g(n)T]  =  S  P  Dg  SH  ,  (A13) 

E[Y5(n)*  ys(n)T]  =  E[v(n)*  y(n)T]  =  al  I  ,  (A14) 

and 

E[v(n)*  v^(n)T]  =  0.  (A  15) 

We  can  now  define  the  so-called  matrix  pencil 

C(X)  =  C-aC.3.  (A  16) 
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where 


C  =  Ryy  -  Cv  I  =  S  P  SH  (A17) 

and 

C5  =  Ryy,  =  SPDg  SH  .  (A18) 

These  correlation  matrices  satisfy  a  so-called  generalized  eigenvalue  problem 

det  (C  -  XCg)  =  0  C  e  =  Cg  e  (A19) 

whose  L  nonzero  eigenvalues  give  the  directions  of  arrival 

Xi  =  e*M  (i  =  1,2,  ...  L)  .  (A 20) 

This  can  be  seen  by  substituting  the  expressions  above  for  C  and  C5  in  the  eigenvalue 
equation.  Premultiplying  both  sides  by  SH  and  removing  the  common  factor  (SHS)P, 
yields 

SH  e  =  X  Dg  SH  e  .  (A21) 

Since  Dg  is  diagonal,  we  have  that  Dg-1  is  equal  to  DgH,  so  that  it  is  easy  to  transfer  Dg  to 
the  other  side  of  the  equation.  If  we  then  define  SHe  as  an  eigenvector  of  Dg,  we  get  (by 
using  the  diagonality  of  Dg)  the  eigenvalue  relation  above,  i.e.,  Xt  =  Dg(i,i). 

In  practice,  we  estimate  the  correlation  matrices  from  the  data,  as  with  MUSIC.  In 
particular,  the  cross-correlation  matrix  is  found  rrom 

Cg  =  Ryyg  =  J-  NI  y(nf  yg(n)T  .  (A22) 

‘ '  n  =  0 

Figures  Dl.l  through  D2  4  are  some  examples  of  results  from  my  ESPRIT  program.  The 
shift  between  the  two  arrays  is  indicated  by  Delta(d),  which  indicates  the  fraction  of  an 
element  spacing  that  separates  the  arrays.  In  practice,  the  array  is  intended  to  consist  of 
matched  pairs  of  elements  rather  than  widely  spaced  arrays.  As  before,  the  element  spacing 
is  X/2.  We  consider  only  the  case  of  perfect  decorrelation,  although  the  method  should 
work  for  perfect  correlation  as  well  (with  suitable  modifications). 

Figures  Dl.l  through  D1.4  compare  MUSIC  and  ESPRIT  in  the  case  of  very  low  noise 
(SNR  =  99  dB).  Both  algorithms  give  exact  results,  suggesting  that  they  both  work 
properly.  (It  should  be  noted,  though,  that  for  the  noiseless  case,  the  generalized 
eigenvalue  algorithm  failed  to  converge  for  certain  values  of  delta,  such  as  0.2,  0.5,  and 
1.0.  This  is  probably  due  to  some  kind  of  singularity  or  ill-posedness  in  the  matrices. 
This  did  not  occur  in  the  noisy  case.)  Figures  D2.1  through  D2.4  present  a  comparison  of 
the  two  in  30  dB  noise  for  the  same  array  separations,  which  seems  to  indicate  that 
ESPRIT  works  even  better  than  does  MUSIC. 
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Appendix  B 

SPATIALLY  AVERAGED  MUSIC 


Spatially  averaged  MUSIC  is  a  variation  of  MUSIC  that  works  in  the  case  of  perfeci 
correlation  between  the  signals.  Essentially,  it  divides  the  array  into  subarrays,  calculates  a 
correlation  matrix  for  each  subarray  in  the  usual  way,  and  then  averages  these  together.  It 
can  be  shown  through  straightforward  substitution  and  algebraic  manipulation  that  the 
singularity  in  the  correlation  matrix  is  removed  by  this  averaging  process*  .  The  different 
subarrays  provide  a  slightly  different  "perspective"  on  the  incoming  signals,  and  this  extra 
degree  of  freedom  overcomes  the  lack  of  full  rank  in  the  unaveraged  correlation  matrix. 
The  price  we  must  pay  is  that  we  are  using  smaller  arrays,  and  it  is  not  immediately  clear 
that  the  decrease  in  resolution  is  acceptable.  This  appendix  presents  some  results  which 
show  that  we  get  performance  in  the  correlated  case  that  rivals  that  of  the  Prony  method 
presented  in  the  main  body  of  this  report. 

Table  B1  presents  a  comparison  of  spatially  averaged  MUSIC  and  Prony.  There  are 
five  elements  in  each  subarray  of  the  spatially  averaged  case.  All  other  quantities  are 
equivalent  to  the  analysis  presented  in  the  report.  "N"  refers  to  the  number  of  snapshots. 
The  numbers  inside  the  table  are  the  angular  thresholds  for  resolution,  i.e.,  the  angular 
separation  in  degrees  at  which  the  two  perfectly  correlated  signals  were  first  resolved.  The 
numbers  above  the  diagonals  are  for  the  averaged  MUSIC  case  and  below  the  diagonals  are 
for  the  Prony  case  as  presented  in  this  report.  Clearly,  the  performance  is  comparable. 


Table  B1 — Comparison  of  Spatially  Averaged  MUSIC  and  Prony 


* 


S.J.  Orphanidis, 


Optimum  Signal  Processing,  Macmillan.  New  York,  N.Y., 
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Appendix  C 

SOME  RECENT  RESEARCH  IN  MUSIC 


The  following  remarks  survey  some  of  the  recent  literature  on  research  into  the 
MUSIC  algorithm. 

MUSIC  is  highly  sensitive  to  system  errors.  Excellent  calibration  is  difficult  to 
accomplish,  and  even  when  possible,  the  system  parameters  will  drift  with  time.  Reference 
Cl  attempts  to  quantify  this  sensitivity  by  a  first-order  analysis,  assuming  a  perfectly 
known  covariance  matrix.  The  basic  conclusion  is  that  increasing  array  aperture  will 
reduce  the  sensitivity,  so  that  there  is  a  tradeoff  between  the  costs  of  calibration  and  of 
increased  aperture. 

A  beamforming  preprocessor  may  increase  the  performance  of  MUSIC  and  of  other 
high-resolution  algorithms.  Benefits  include  reduced  computation,  improved  performance 
in  the  presence  of  spatially  colored  noise,  and  increased  resolution.  Reference  C2  provides 
an  in-depth,  quantitative  analysis  for  closely  spaced,  uncorrelated  targets. 

Root-MUSIC  is  a  variation  of  MUSIC  that  has  superior  performance  in  the  case  of  a 
linear,  equispaced  array.  This  method  determines  the  direction-of-arrivals  by  solving  a 
polynomial  formed  from  the  noise  subspace,  rather  than  scanning  numerically  over  all 
DOAs  to  find  the  orthogonality  between  signal  and  noise  subspaces  described  in  this 
report.  Reference  C3  explains  and  quantifies  Root-MUSIC  with  a  perturbation  analysis  on 
the  roots  of  the  polynomial. 

In  Reference  C4  an  asymptotic  analysis  is  presented  for  a  generalized  version  of  the 
smoothing  scheme  for  coherent  sources  discussed  in  Appendix  B  of  this  report.  This 
generalization  is  called  forward-backward,  smoothing  as  compared  to  the  forward  only 
smoothing  in  Appendix  B.  Forward  only  smoothing  requires  2K  elements  to  detect  K 
coherent  DOAs,  whereas  forward-backward  smoothing  requires  only  [3K/2]  elements.  In 
Pillai's  analysis,  a  uniform  array  is  assumed,  which  is  a  restriction  on  the  arbitrary  array 
geometry  of  which  MUSIC  is  capable.  A  key  result  is  that  about  (1/Mcod)2  *  1  more 
snapshots  are  required  to  resolve  two  equipowered  closely  spaced  coherent  plane  waves 
than  are  required  for  uncorrelated  waves,  where  M  is  the  number  of  elements  and  0)<j  is  the 
angular  separation.  Statistics  for  the  unsmoothed  and  forward-only  smoothed  cases  result 
as  special  cases  of  Pillai's  general  analysis. 

Reference  C5  provides  analytic  results  for  the  asymptotic  efficiency  of  MUSIC  at  high 
SNR  with  multiple,  uncorrelated  targets.  That  is,  when  the  number  of  snapshots 
approaches  infinity,  the  variances  of  the  DOAs  approach  the  Cramer-Rao  lower  bound. 

Minimum-Norm  is  another  leading  contender  besides  MUSIC  for  array  high  resolution. 
An  asymptotic  statistical  analysis  of  these  two  algorithms  is  provided  in  Reference  C6. 
Results  suggest  that  Min-Norm  has  a  resolution  threshold  at  a  lower  SNR  than  MUSIC. 

Reference  C7  presents  an  analysis  of  the  spatially  averaged  method  for  correlated  signals 
discussed  in  Appendix  B. 
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